/* -----------------------------------------------------------------------------
 * dicom.cpp 
 * -----------------------------------------------------------------------------
 * Copyright (c) 2015-2017 Blaine Rister et al., see LICENSE for details.
 * -----------------------------------------------------------------------------
 * C-language wrapper for the DCMTK library.
 * -----------------------------------------------------------------------------
 */

/* SIFT3D includes */
#include "imutil.h"
#include "immacros.h"
#include "dicom.h"

#ifndef SIFT3D_WITH_DICOM
/* Return error messages if this was not compiled with DICOM support. */ 

static int dcm_error_message() {
        SIFT3D_ERR("dcm_error_message: SIFT3D was not compiled with DICOM "
                "support!\n");
        return SIFT3D_WRAPPER_NOT_COMPILED;
}

int read_dcm(const char *path, Image *const im) {
        return dcm_error_message();
}

int read_dcm_dir(const char *path, Image *const im) {
        return dcm_error_message();
}

int write_dcm(const char *path, const Image *const im, 
        const Dcm_meta *const meta, const float max_val) {
        return dcm_error_message();
}

int write_dcm_dir(const char *path, const Image *const im, 
        const Dcm_meta *const meta) {
        return dcm_error_message();
}

#else

/*----------------Include the very picky DCMTK----------------*/
#include "dcmtk/config/osconfig.h"    /* make sure OS specific configuration is included first */

#define INCLUDE_CSTDIO
#define INCLUDE_CSTRING
#include "dcmtk/ofstd/ofstdinc.h"

#ifdef HAVE_GUSI_H
#include <GUSI.h>
#endif

#include "dcmtk/dcmdata/dctk.h"          /* for various dcmdata headers */
#include "dcmtk/dcmdata/cmdlnarg.h"      /* for prepareCmdLineArgs */
#include "dcmtk/dcmdata/dcuid.h"         /* for dcmtk version name */

#include "dcmtk/ofstd/ofconapp.h"        /* for OFConsoleApplication */
#include "dcmtk/ofstd/ofcmdln.h"         /* for OFCommandLine */

#include "dcmtk/oflog/oflog.h"           /* for OFLogger */

#include "dcmtk/dcmimgle/dcmimage.h"     /* for DicomImage */
#include "dcmtk/dcmimgle/diutils.h"     /* for DIPixel */
#include "dcmtk/dcmimage/diregist.h"     /* include to support color images */
#include "dcmtk/dcmdata/dcrledrg.h"      /* for DcmRLEDecoderRegistration */

#include "dcmtk/dcmjpeg/djdecode.h"      /* for dcmjpeg decoders */
#include "dcmtk/dcmjpeg/dipijpeg.h"      /* for dcmimage JPEG plugin */

#include "dcmtk/dcmjpeg/djencode.h" /* for JPEG encoding */
#include "dcmtk/dcmjpeg/djrplol.h"  /* for DJ_RPLossless */

#include "dcmtk/dcmseg/segment.h" /* Dicom segmentations */

#ifdef WITH_ZLIB
#include <zlib.h>          /* for zlibVersion() */
#endif
/*---------------------------------------------------------*/

/* Other includes */
#include <algorithm>
#include <memory>
#include <vector>
#include <cmath>
#include <cfloat>
#include <cstdlib>
#include <stdint.h>
#include <dirent.h>

/* Macro to call a C++ function and catch any exceptions it throws,
 * returning SIFT3D_FAILURE when an exception is caught. The return value is
 * stored in ret. */
#define CATCH_EXCEPTIONS(ret, tag, fun, ...) \
        try { \
                ret = (fun)( __VA_ARGS__ ); \
        } catch (std::exception &e) { \
                SIFT3D_ERR("%s: %s\n", tag, e.what()); \
                ret = SIFT3D_FAILURE; \
        } catch (...) { \
                SIFT3D_ERR("%s: unexpected exception \n", tag); \
                ret = SIFT3D_FAILURE; \
        } \

/* Format strings for DICOM tags */
static const char *pixelSpacingFmt = "%lf\\%lf"; // Pixel Spacing
static const char *imPosPatientFmt = "%f\\%f\\%f"; // ImagePositionPatient
static const char *imOriPatientFmt = "%f\\%f\\%f\\%f\\%f\\%f"; 
        // ImageOrientationPatient

/* File separator in string form */
const std::string sepStr(1, SIFT3D_FILE_SEP);

/* Dicom parameters */
static const unsigned int dcm_bit_width = 8; // Bits per pixel

/* DICOM metadata defaults */
static const char *default_patient_name = "DefaultSIFT3DPatient";
static const char *default_series_descrip = "Series generated by SIFT3D";
static const char *default_patient_id = "DefaultSIFT3DPatientID";
static const char default_instance_num = 1;

/* Types of DICOMs which necessitate special processing */
enum image_type {
        DSO,
        PET,
        OTHER
};

/* Helper declarations */
class Dicom;
static bool isLittleEndian(void);
static void default_Dcm_meta(Dcm_meta *const meta);
static int load_file(const char *path, DcmFileFormat &fileFormat);
static int read_dcm_cpp(const char *path, Image *const im);
static int read_dcm_img(const Dicom &dicom, Image *const im);
static int read_dso(const char *imDir, Dicom &dso, 
                            Image *const mask);
static int cvec_max_abs(const Cvec *v, float *const val, int *const pos);
static int read_dcm_dir_meta(const char *path, std::vector<Dicom> &dicoms);
static int read_dcm_dir_cpp(const char *path, Image *const im);
static int write_dcm_cpp(const char *path, const Image *const im,
        const Dcm_meta *const meta, const float max_val);
static int write_dcm_dir_cpp(const char *path, const Image *const im,
        const Dcm_meta *const meta);
static void set_meta_defaults(const Dcm_meta *const meta, 
        Dcm_meta *const meta_new);
static int dcm_resize_im(const std::vector<Dicom> &dicoms, Image *const im);
static int write_subvolume(const Dicom &dicom, Image *const main, 
        const int mainOffset, Image *const sub, const int start, 
        const int end);

/* Helper class to store DICOM data. */
class Dicom {
private:
        int axes[2]; // Data axes, e.g. {0, 1} means x, y
        int axisSigns[2]; // Axis signs, negative means go backwards
        std::string filename; // DICOM file name
        std::string classUID;
        std::string seriesUID; // Series UID 
        std::string instanceUID; // Instance UID
        double sortCoord; // Coordinate by which the series is sorted. Usually z
        double ux, uy, uz; // Voxel spacing in real-world coordinates
        int nx, ny, nz, nc; // Image dimensions
        int sortAxis; // The axis corresponding to sortCoord
        bool valid; // Data validity 

public:

        /* Data is initially invalid */
        Dicom() : valid(false) {};

        ~Dicom() {};

        /* Load a file */
        Dicom(const char *filename, const int defaults3D=0);

        /* Get the multiplier to convert a PET image to SUV */
        double PET_SUV_multiplier(void) const;

        /* Get a dimension by index */
        int getDim(const int idx) const {
                switch (idx) {
                case 0:
                        return getNx();
                case 1:
                        return getNy();
                case 2:
                        return getNz();
                }

                return -1;
        }

        /* Get a unit by index */
        int getUnit(const int idx) const {
                switch (idx) {
                case 0:
                        return getUx();
                case 1:
                        return getUy();
                case 2:
                        return getUz();
                }

                return -1;
        }

        /* Get the image axis by index */
        int getAxes(const int idx) const {
                return idx < 0 || idx > 2 ? -1 : axes[idx];
        }

        /* Get the axis sign by index */
        int getAxisSign(const int idx) const {
                return idx < 0 || idx > 2 ? -1 : axisSigns[idx];
        }

        /* Get the x-dimension */
        int getNx(void) const {
                return nx;
        }

        /* Get the y-dimension */
        int getNy(void) const {
                return ny;
        }

        /* Get the z-dimension */
        int getNz(void) const {
                return nz;
        }

        /* Get the number of channels */
        int getNc(void) const {
                return nc;
        } 

        /* Get the sorting axis */
        int getSortAxis(void) const {
                return sortAxis;
        }

        /* Get the sorting coordinate (usually z) */
        double getSortCoord(void) const {
                return sortCoord;
        }

        /* Get the units in the sorting dimension */
        int getSortUnit(void) const {
                return getUnit(getSortAxis());
        }

        /* Get the x-spacing */
        double getUx(void) const {
                return ux;
        }

        /* Get the y-spacing */
        double getUy(void) const {
                return uy;
        }

        /* Get the z-spacing */
        double getUz(void) const {
                return uz;
        }

        /* Check whether or not the data is valid */
        bool isValid(void) const {
                return valid;
        }

        /* Get the file name */
        const char *name(void) const {
                return filename.c_str();
        }

        /* Sort by z position */
        bool operator < (const Dicom &dicom) const {
                return getSortCoord() < dicom.getSortCoord();
        }

        /* Check if another DICOM file is from the same series */
        bool eqSeries(const Dicom &dicom) const {
                return !seriesUID.compare(dicom.seriesUID);
        }

        /* Check for a matching SOPInstanceUID */
        bool eqInstance(const char *uid) const {
                return !instanceUID.compare(uid);
        }

        /* Return an image type, based on the classUID */
        enum image_type getType() const {
                if (!classUID.compare(UID_SegmentationStorage)) {
                        return DSO;
                } else if (
                        !classUID.compare(
                                UID_PositronEmissionTomographyImageStorage) ||
                        !classUID.compare(
                                UID_LegacyConvertedEnhancedPETImageStorage) ||
                        !classUID.compare(UID_EnhancedPETImageStorage)) {
                        return PET;
                } else {
                        return OTHER;
                }
        }
};

/* Read a DICOM file with DCMTK. */
static int load_file(const char *path, DcmFileFormat &fileFormat) {

        // Load the image as a DcmFileFormat 
        OFCondition status = fileFormat.loadFile(path);
        if (status.bad()) {
                SIFT3D_ERR("load_file: failed to read DICOM file %s (%s)\n",
                        path, status.text());
                return SIFT3D_FAILURE;
        }

        return SIFT3D_SUCCESS;
}

/* Get the single maximum component from a Cvec, by absolue value. Returns the 
 * element and its position. Returns SIFT3D_SUCCESS if the maximum is unique,
 * SIFT3D_FAILURE otherwise. */
static int cvec_max_abs(const Cvec *v, float *const val, int *const pos) {

        int k;
        const float vals[] = {v->x, v->y, v->z};

        // Find the maximum
        float maxDiff = 0.f;
        float maxAbsVal = -1;
        for (k = 0; k < 3; k++) {

                const float thisVal = vals[k];
                const float thisAbsVal = fabsf(thisVal);
        
                // Check for maximum. If so, write into val and pos
                if (thisAbsVal < maxAbsVal)
                        continue;
                maxDiff = thisAbsVal - maxAbsVal;
                maxAbsVal = thisAbsVal;
                *val = thisVal;
                *pos = k;
        }

        // Test for unique maxima
        return maxDiff > 1e-2 ? SIFT3D_SUCCESS : SIFT3D_FAILURE;
}

/* Load the data from a DICOM file. If defaults3D is true, substitute default
 * values for 3D positioning information. Use this to load a single file even
 * if key metadata is missing, e.g. for 2D images. */
Dicom::Dicom(const char *path, const int defaults3D) : filename(path), 
        valid(false) {

        // Read the file
        DcmFileFormat fileFormat;
        if (load_file(path, fileFormat))
                return;
        DcmDataset *const data = fileFormat.getDataset();

        // Get the SOPClass UID
        const char *classUIDStr;
        OFCondition status = data->findAndGetString(DCM_SOPClassUID, 
                                                    classUIDStr);
        if (status.bad() || classUIDStr == NULL) {
                SIFT3D_ERR("Dicom.Dicom: failed to get SOPClassUID "
                        "from file %s (%s)\n", path, status.text());
                return;
        }
        classUID = std::string(classUIDStr); 

        // Get the series UID 
        const char *seriesUIDStr;
        status = data->findAndGetString(DCM_SeriesInstanceUID, seriesUIDStr);
        if (status.bad() || seriesUIDStr == NULL) {
                SIFT3D_ERR("Dicom.Dicom: failed to get SeriesInstanceUID "
                        "from file %s (%s)\n", path, status.text());
                return;
        }
        seriesUID = std::string(seriesUIDStr); 

        // Get the instance UID 
        const char *instanceUIDStr;
        status = data->findAndGetString(DCM_SOPInstanceUID, instanceUIDStr);
        if (status.bad() || instanceUIDStr == NULL) {
                SIFT3D_ERR("Dicom.Dicom: failed to get SOPInstanceUID "
                        "from file %s (%s)\n", path, status.text());
                return;
        }
        instanceUID = std::string(instanceUIDStr); 

        // Get the z coordinate
        if (getType() == DSO) {
                // DSOs don't always have coordinates
                sortCoord = -1;

                // DSOs don't always have units
                ux = uy = uz = -1.0;
        } else {
#if 0
                // Read the patient position
                const char *patientPosStr;
                status = data->findAndGetString(DCM_PatientPosition, 
                                                patientPosStr);
                if (status.bad() || patientPosStr == NULL) {
                        SIFT3D_ERR("Dicom.Dicom: failed to get " 
                                "PatientPosition from file %s (%s)\n", path, 
                                status.text());
                        return;
                }

                // Interpret the patient position to give the sign of the z axis
                double zSign;
                switch (patientPosStr[0]) {
                case 'H':
                        zSign = -1.0;
                        break;
                case 'F':
                        zSign = 1.0;
                        break;
                default:
                        SIFT3D_ERR("Dicom.Dicom: unrecognized patient "
                                "position: %d\n", patientPosStr);
                        return;
                }
#else
                //TODO: Is this needed?
                const double zSign = 1.0;
#endif

                // Read the image position patient vector
                const char *imPosPatientStr;
                status = data->findAndGetString(DCM_ImagePositionPatient, 
                        imPosPatientStr);
                if (status.bad() || imPosPatientStr == NULL) {
                        SIFT3D_ERR("Dicom.Dicom: failed to get "
                        "ImagePositionPatient from file %s (%s)\n, defaulting"
                        " to zeros. \n", path, 
                        status.text());
                        if (defaults3D) {
                                imPosPatientStr = "0\\0\\0";
                        } else return;
                }

                // Parse the image position patient vector
                Cvec imPos;
                if (sscanf(imPosPatientStr, imPosPatientFmt, &imPos.x, &imPos.y, 
                        &imPos.z) != 3) {
                        SIFT3D_ERR("Dicom.Dicom: failed to parse "
                                "ImagePositionPatient value %s from file %s\n", 
                                imPosPatientStr, path);
                        if (defaults3D) {
                                imPos = {0, 0, 0};
                        } else return;
                }

                // Read the image orientation patient vector
                const char *imOriPatientStr; 
                status = data->findAndGetString(DCM_ImageOrientationPatient,
                        imOriPatientStr);
                if (status.bad() || imOriPatientStr == NULL) {
                        SIFT3D_ERR("Dicom.Dicom: failed to get "
                        "ImageOrientationPatient from file %s (%s)\n", path, 
                        status.text());
                        return;
                }

                // Parse the image orientation patient vectors
                Cvec imOri1, imOri2;
                if (sscanf(imOriPatientStr, imOriPatientFmt, 
                        &imOri1.x, &imOri1.y, &imOri1.z, &imOri2.x, &imOri2.y,
                        &imOri2.z) != 6) {
                        SIFT3D_ERR("Dicom.Dicom: failed to parse "
                                "ImageOrientationPatient value %s from file %s\n", 
                                imOriPatientStr, path);
                        return;
                }

                // Take the cross-product of orientation vectors to get the 
                // image normal
                Cvec normal; 
                SIFT3D_CVEC_CROSS(&imOri1, &imOri2, &normal);

                // Project the patient position on the image normal, to get the
                // sorting coordinate
                sortCoord = SIFT3D_CVEC_DOT(&imPos, &normal);

                // Get the image axes
                float oriVals[2];
                if (cvec_max_abs(&imOri1, oriVals, axes) || 
                        cvec_max_abs(&imOri2, oriVals + 1, axes + 1)) {
                        SIFT3D_ERR("Dicom.Dicom: unsupported format for "
                                "imageOrientationPatient %s from file %s\n", 
                                imOriPatientStr, path);
                        return;
                }

                // Get the sorting axis
                for (int k = 0; k < 3; k++) {
                        if (axes[0] == k || axes[1] == k)
                                continue;
                        sortAxis = k;
                        break;
                }

                // Compute the signs of the axes
                for (int k = 0; k < 2; k++) {
                        axisSigns[k] = oriVals[k] >= 0 ? 1 : -1;
                }

                // Read the pixel spacing
                double spacing[2];
                const char *pixelSpacingStr;
                status = data->findAndGetString(DCM_PixelSpacing, 
                                                pixelSpacingStr);
                if (status.bad()) {
                        SIFT3D_ERR("Dicom.Dicom: failed to get pixel spacing "
                                "from file %s (%s)\n", path, status.text());
                        return;
                }
                if (sscanf(pixelSpacingStr, pixelSpacingFmt, spacing, 
                        spacing + 1) != 2) {
                        SIFT3D_ERR("Dicom.Dicom: unable to parse pixel "
                                "spacing from file %s \n", path);
                        return;
                }
                if (spacing[0] <= 0.0 || spacing[1] <= 0.0) {
                        SIFT3D_ERR("Dicom.Dicom: file %s has invalid pixel "
                                "spacing [%f, %f]\n", path, ux, uy);
                        return;
                }

                // Convert the pixel spacing into units based on the image axes
                double *units = &ux;
                for (int k = 0; k < 2; k++) {
                       units[axes[k]] = spacing[k]; 
                }

                // Read the slice thickness 
                Float64 sliceThickness;
                status = data->findAndGetFloat64(DCM_SliceThickness, 
                                                 sliceThickness);
                if (!status.good()) {
                        SIFT3D_ERR("Dicom.Dicom: failed to get slice "
                                "thickness from file %s (%s)\n", path, 
                                status.text());
                        return;
                }
                if (sliceThickness <= 0.0) {
                        SIFT3D_ERR("Dicom.Dicom: file %s has invalid slice "
                                "thickness: %f \n", path, uz);
                        return;
                }

                // Use the slice thickness as the units for the perpendicular 
                // axis
                units[sortAxis] = sliceThickness;
        }

        // Load the DicomImage object
        DicomImage dicomImage(path);
        if (dicomImage.getStatus() != EIS_Normal) {
                SIFT3D_ERR("Dicom.Dicom: failed to open image %s (%s)\n",
                        path, 
                        DicomImage::getString(dicomImage.getStatus()));
                return;
        }

        // Check for color images
        if (!dicomImage.isMonochrome()) {
                SIFT3D_ERR("Dicom.Dicom: reading of color DICOM images is "
                        "not supported at this time \n");
                return;
        }
        nc = 1;

        // Read the dimensions
        int *dims = &nx;
        dims[axes[0]] = dicomImage.getWidth();
        dims[axes[1]] = dicomImage.getHeight();
        dims[sortAxis] = dicomImage.getFrameCount();
        if (nx < 1 || ny < 1 || nz < 1) {
                SIFT3D_ERR("Dicom.Dicom: invalid dimensions for file %s "
                        "(%d, %d, %d)\n", path, nx, ny, nz);
                return;
        }

        // Set the window 
        dicomImage.setMinMaxWindow();

        valid = true;
}

/* Check the endianness of the machine. Returns true if the machine is little-
 * endian, false otherwise. */
static bool isLittleEndian(void) {
        volatile uint16_t i = 0x0123;
        return ((uint8_t *) &i)[0] == 0x23;
}

/* Set a Dcm_meta struct to default values. Generates new UIDs. */
static void default_Dcm_meta(Dcm_meta *const meta) {
        meta->patient_name = default_patient_name;
        meta->patient_id = default_patient_id;
        meta->series_descrip = default_series_descrip;
        dcmGenerateUniqueIdentifier(meta->study_uid, SITE_STUDY_UID_ROOT);
        dcmGenerateUniqueIdentifier(meta->series_uid, SITE_SERIES_UID_ROOT);
        dcmGenerateUniqueIdentifier(meta->instance_uid, SITE_INSTANCE_UID_ROOT); 
        meta->instance_num = default_instance_num;
}

/* Parse a time (TM) string from a Dicom file. */ 
int parseTM(const char *const tm, double *const time) {

        // Buffers to hold the sub-strings
        char hh[3];
        char mm[3];

        // Subdivide into the TM into hours, minutes and seconds
        memcpy(hh, tm, 2);
        memcpy(mm, tm + 2, 2);
        hh[2] = mm[2] = '\0';
        const char *const ss = tm + 4;

        // Parse the strings
        errno = 0;
        const int hours = atoi(hh);
        const int minutes = atoi(mm);
        const double seconds = strtod(ss, NULL);
        if (errno) return SIFT3D_FAILURE;

        // Add up the totals
        *time = hours * 60.0 * 60.0 + minutes * 60.0 + seconds;

        return SIFT3D_SUCCESS;
}

/* Get the multiplier to convert a PET scan to SUV values. Returns negative if
 * an error has occured. */
double Dicom::PET_SUV_multiplier(void) const {

        // Read the file
        const char *path = filename.c_str();
        DcmFileFormat fileFormat;
        if (load_file(path, fileFormat))
                return -1.0;
        DcmDataset *const data = fileFormat.getDataset();

        // Read the body weight (0010, 1010)
        Float64 weight;
        OFCondition status = data->findAndGetFloat64(DCM_PatientWeight, 
                weight);
        if (status.bad()) {
                SIFT3D_ERR("Dicom.PET_SUV_multiplier: failed to get "
                        "PatientWeight (0010, 1010) from file %s (%s)\n", 
                        path, status.text());
                return -1.0;
        }

        // Read the radiopharmeceutical dosage (0018, 1074)
        Float64 dosage;
        status = data->findAndGetFloat64(DCM_RadionuclideTotalDose,
                dosage, 0, OFTrue);
        if (status.bad()) {
                SIFT3D_ERR("Dicom.PET_SUV_multiplier: failed to get "
                        "RadionuclideTotalDose (0018, 1074) from file %s "
                        "(%s)\n", path, status.text());
                return -1.0;
        }

        // Read the radiopharmeceutical injection time (0018, 1072)
        // DCM_RadiopharmaceuticalStartTime (TM)
        const char *timeStr;
        status = data->findAndGetString(DCM_RadiopharmaceuticalStartTime, 
                timeStr, OFTrue);
        if (status.bad() || timeStr == NULL) {
                SIFT3D_ERR("Dicom.PET_SUV_multiplier: failed to get "
                        "RadiopharmeceuticalStartTime (0018, 1072) from file %s"
                        " (%s)\n", path, status.text());
                return -1.0;
        }

        // Parse the injection time
        Float64 iv_time;
        if (parseTM(timeStr, &iv_time)) {
                SIFT3D_ERR("Dicom.PET_SUV_multipler: failed to parse "
                "RadiopharmeceuticalStartTime (0018, 1072) value %s from file "
                "%s\n", timeStr, path);
                return -1.0;
        }

        // Read the image acquistion time (0008, 0032)
        status = data->findAndGetString(DCM_AcquisitionTime, timeStr);
        if (status.bad() || timeStr == NULL) {
                SIFT3D_ERR("Dicom.PET_SUV_multiplier: failed to get "
                        "AcquisitionTime (0008, 0032) from file %s"
                        " (%s)\n", path, status.text());
                return -1.0;
        }

        // Parse the image time
        Float64 image_time;
        if (parseTM(timeStr, &image_time)) {
                SIFT3D_ERR("Dicom.PET_SUV_multipler: filed to parse "
                        "AcquisitionTime (0008, 0032) value %s from file %s\n", 
                        timeStr, path);
                return -1.0;
        }

        // Read the radiopharmeceutical half-life (0018, 1075)
        Float64 half_life;
        status = data->findAndGetFloat64(DCM_RadionuclideHalfLife, half_life, 0,
                OFTrue);
        if (status.bad()) {
                SIFT3D_ERR("Dicom.PET_SUV_multiplier: failed to get "
                        "RadionuclideHalfLife (0018, 1075) from file %s (%s)\n",
                        path, status.text());
                return -1.0;
        }

        // Compute the elapsed time
        double elapsed_time = iv_time - image_time; 
        if (elapsed_time < 0) {
                const double day_seconds = 24.0 * 60.0 * 60.0;
                elapsed_time = elapsed_time + day_seconds;
        }

        // Compute the time-adjusted dosage
        const double adjusted_dosage = dosage * 
                pow(2.0, -elapsed_time / half_life);

        // Compute the SUV multiplier
        return weight / adjusted_dosage;
}

/* Read a DICOM file into an Image struct. If the file is a DSO, the 
 * referenced DICOM image must be in the same directory. */
int read_dcm(const char *path, Image *const im) {

        int ret;

        CATCH_EXCEPTIONS(ret, "read_dcm", read_dcm_cpp, path, im);

        return ret;
}

/* Read all of the DICOM files from a directory into an Image struct. Slices 
 * must be ordered alphanumerically, starting with z = 0. */
int read_dcm_dir(const char *path, Image *const im) {

        int ret;

        CATCH_EXCEPTIONS(ret, "read_dcm_dir", read_dcm_dir_cpp, path, im);

        return ret;
}

/* Write an Image struct into a DICOM file. 
 * Inputs: 
 *      path - File name
 *      im - Image data
 *      meta - Dicom metadata (or NULL for default values)
 *      max_val - The maximum value of the image, used for scaling. If set
 *              to a negative number, this functions computes the maximum value
 *              from this image.
 *
 * Returns SIFT3D_SUCCESS on success, SIFT3D_FAILURE otherwise.
 */
int write_dcm(const char *path, const Image *const im, 
        const Dcm_meta *const meta, const float max_val) {

        int ret;

        CATCH_EXCEPTIONS(ret, "write_dcm", write_dcm_cpp, path, im, meta, 
                max_val);

        return ret;
}

/* Write an Image struct into a directory of DICOM files.
 * Inputs: 
 *      path - File name
 *      im - Image data
 *      meta - Dicom metadata (or NULL for default values)
 *
 * Returns SIFT3D_SUCCESS on success, SIFT3D_FAILURE otherwise.
 */
int write_dcm_dir(const char *path, const Image *const im, 
        const Dcm_meta *const meta) {

        int ret;

        CATCH_EXCEPTIONS(ret, "write_dcm_dir", write_dcm_dir_cpp, path, im, 
                meta);

        return ret;
}

/* Helper function to read a DICOM file using C++ */
static int read_dcm_cpp(const char *path, Image *const im) {

        // Read the image metadata
	Dicom dicom(path, 1);
        if (!dicom.isValid())
                return SIFT3D_FAILURE;

        // Check if this is a segmentation or normal image
        if (dicom.getType() == DSO) {
                // Look for images only in the same directory
                char *pathDir = im_get_parent_dir(path);
                int ret = read_dso(pathDir, dicom, im);
                free(pathDir);
                return ret;
         }
       
        // Otherwise directly read the image data 
        return read_dcm_img(dicom, im);
}

/* Helper function to read DICOM image data */
static int read_dcm_img(const Dicom &dicom, Image *const im) {

        int offsets[] = {0, 0, 0};
        int signs[] = {1, 1, 1};
        const void *data;
        const DiMonoPixel *pixels;
	uint32_t shift;
	int x, y, z, depth, k;
        const int bufNBits = 32;

        // Image data strides
        const int y_stride = dicom.getNx(); 
        const int z_stride = dicom.getNx() * dicom.getNy();

	// Initialize JPEG decoders
	DJDecoderRegistration::registerCodecs();

        // Initialize the DicomImage object
        const char *path = dicom.name();
	DicomImage dicomImage(path);
        if (dicomImage.getStatus() != EIS_Normal) {
                SIFT3D_ERR("read_dcm_cpp: failed to open image %s (%s)\n",
                        path, DicomImage::getString(dicomImage.getStatus()));
		goto read_dcm_img_quit;
        }

        // Initialize the image fields
        im->nx = dicom.getNx();
        im->ny = dicom.getNy();
        im->nz = dicom.getNz();
        im->nc = dicom.getNc();
        im->ux = dicom.getUx();
        im->uy = dicom.getUy();
        im->uz = dicom.getUz();

        // Resize the output
        im_default_stride(im);
        if (im_resize(im))
		goto read_dcm_img_quit;

        // Format the offsets and signs for the main volume
        for (k = 0; k < 2; k++) {
                if (dicom.getAxisSign(k) > 0) 
                        continue;
                const int axisIdx = dicom.getAxes(k);
                signs[axisIdx] = -1;
                offsets[axisIdx] = dicom.getDim(axisIdx) - 1;
        }

        // Get the vendor-independent intermediate pixel data
        pixels = (const DiMonoPixel *) dicomImage.getInterData();
        if (pixels == NULL) {
                SIFT3D_ERR("read_dcm_img: failed to get intermediate data for "
                        "%s\n", path);
                goto read_dcm_img_quit;
        }
        depth = pixels->getBits();

        // Macro to copy the data
#define COPY_DATA(type) { \
        SIFT3D_IM_LOOP_START(im, x, y, z) \
                SIFT3D_IM_GET_VOX(im, x * signs[0] + offsets[0], \
                        y * signs[1] + offsets[1], \
                        z * signs[2] + offsets[2], \
                        0) = (float) *((type *) data + x + y * y_stride + \
                                z * z_stride);\
        SIFT3D_IM_LOOP_END \
        }\

        // Choose the appropriate data type and copy the data
        data = pixels->getData(); 
        switch (pixels->getRepresentation()) {
        case EPR_Uint8:
                COPY_DATA(uint8_t)
                break;
        case EPR_Uint16:
                COPY_DATA(uint16_t)
                break;
        case EPR_Uint32:
                COPY_DATA(uint32_t)
                break;
        case EPR_Sint8:
                COPY_DATA(int8_t)
                break;
        case EPR_Sint16:
                COPY_DATA(int16_t)
                break;
        case EPR_Sint32:
                COPY_DATA(int32_t)
                break;
        default:
                SIFT3D_ERR("read_dcm_img: unrecognized pixel representation "
                        "for %s\n", path);
                goto read_dcm_img_quit;
        }
#undef COPY_DATA

#if 0
        // Get the bit depth of the image
        depth = dicomImage.Depth();
        if (depth > bufNBits) 
                SIFT3D_ERR("read_dcm_cpp: buffer is insufficiently wide "
                        "for %d-bit data of image %s \n", depth, path);
		goto read_dcm_img_quit;
        }

        // Get the number of bits by which we need to shift the 32-bit data, to
        // recover the original resolution (DICOM uses Big-endian encoding)
        shift = isLittleEndian() ? static_cast<uint32_t>(bufNBits - depth) : 0;

        // Read each frame
        for (int i = 0; i < im->nz; i++) { 

                int x, y, z;

                const int x_start = 0;
                const int y_start = 0;
                const int z_start = i;
                const int x_end = im->nx - 1;
                const int y_end = im->ny - 1;
                const int z_end = z_start;

                // Get a pointer to the data, rendered as a 32-bit int
                const uint32_t *const frameData = 
                        static_cast<const uint32_t *const>(
                                dicomImage.getOutputData(
                                        static_cast<int>(bufNBits), i));
                if (frameData == NULL) {
                        SIFT3D_ERR("read_dcm_cpp: could not get data from "
                                "image %s frame %d (%s)\n", path, i, 
                                DicomImage::getString(dicomImage.getStatus()));
			goto read_dcm_img_quit;
                }

                // Copy the frame
                SIFT3D_IM_LOOP_LIMITED_START(im, x, y, z, x_start, x_end,
                        y_start, y_end, z_start, z_end)

                        // Get the voxel and shift it to match the original
                        // magnitude 
                        const uint32_t vox =
                                frameData[x + y * im->nx] >> shift;

                        // Convert to float and write to the output image
                        SIFT3D_IM_GET_VOX(im, x, y, z, 0) =
                                static_cast<float>(vox);

                SIFT3D_IM_LOOP_END
        }
#endif

	// Clean up
	DJDecoderRegistration::cleanup();

        // Modality-specific post-processing
        if (dicom.getType() == PET) {

                const double suv_factor = dicom.PET_SUV_multiplier();

                // Check for errors in SUV computation
                if (suv_factor < 0)
                        return SIFT3D_FAILURE;

                // Scale all the values
                SIFT3D_IM_LOOP_START(im, x, y, z)
                        SIFT3D_IM_GET_VOX(im, x, y, z, 0) *= suv_factor;
                SIFT3D_IM_LOOP_END
        }

        return SIFT3D_SUCCESS;

read_dcm_img_quit:
	DJDecoderRegistration::cleanup();
	return SIFT3D_FAILURE;
}

/* Read a DICOM Segmentation Object (DSO) mask. 
 *
 * Parameters:
 *      imDir - The path of the directory containing DICOM image files.
 *      dsoPath - The path of the DSO file.
 *      mask - The image in which to write a binary mask.
 *
 * Returns SIFT3D_SUCCESS on success, SIFT3D_FAILURE otherwise.
 */
static int read_dso(const char *imDir, Dicom &dso, 
                            Image *const mask) {

        // Read the DSO
        DcmFileFormat fileFormat;
        const char *dsoPath = dso.name();
        if (load_file(dso.name(), fileFormat))
                return SIFT3D_FAILURE;
        DcmDataset *const dso_data = fileFormat.getDataset();

        // Initialize the DcmSegmentation object
        OFCondition status;
        OFunique_ptr<DcmSegmentation> dcmSegmentation;
        {
                DcmSegmentation *p;
                status = DcmSegmentation::loadDataset(*dso_data, p);
                dcmSegmentation.reset(p);
        }
        if (status.bad()) {
                SIFT3D_ERR("read_dso: failed to initialize the "
                        "DcmSegmentation object for file %s (%s)\n", dsoPath,
                        status.text());
                return SIFT3D_FAILURE;
        }

        // Ensure we have only one segment
        const int num_segments = dcmSegmentation->getNumberOfSegments();
        if (num_segments != 1) {
                SIFT3D_ERR("read_dso: unsupported number of segments "
                        " in file %s: %d\n", dsoPath, num_segments);
                return SIFT3D_FAILURE;
        }

        // Check how many frames we have
        const int nz = dcmSegmentation->getNumberOfFrames();

        // Get the referenced series items
        OFVector<IODSeriesAndInstanceReferenceMacro::ReferencedSeriesItem*>&
                series_items = dcmSegmentation->getCommonInstanceReference()
                        .getReferencedSeriesItems();

        // Ensure there is only one series
        const int num_series = series_items.size();
        if (num_series != 1) {
                SIFT3D_ERR("read_dso: unsupported number of referenced "
                        "series: %d \n", num_series);
                return SIFT3D_FAILURE;
        }

        // Get the referenced instance items
        OFVector<SOPInstanceReferenceMacro*> &ref_instances =
                series_items[0]->getReferencedInstanceItems();

        // Read the image metadata
        std::vector<Dicom> images;
        if (read_dcm_dir_meta(imDir, images))
                return SIFT3D_FAILURE;

        // Read the mask DICOM image data
        Image maskData;
        init_im(&maskData);
        if (read_dcm_img(dso, &maskData)) {
                im_free(&maskData);
                return SIFT3D_FAILURE;
        }

        // Ensure we have a single frame for each instance
        if (ref_instances.size() != maskData.nz) {
                SIFT3D_ERR("read_dso: DSO has %lu frames but %d segments \n",
                        ref_instances.size(), maskData.nz);
                im_free(&maskData);
                return SIFT3D_FAILURE;
        }

        // Initialize and zero the mask
        if (dcm_resize_im(images, mask)) {
                im_free(&maskData);
                return SIFT3D_FAILURE;
        }
        im_zero(mask);

        // Ensure the x,y dimensions match
        if (mask->nx != maskData.nx || mask->ny != maskData.ny) {
                SIFT3D_ERR("read_dso: Mask has frame dimensions [%d x %d] "
                        "while image volume has dimensions [%d x %d] \n",
                        maskData.nx, maskData.ny, mask->nx, mask->ny);
                im_free(&maskData);
                return SIFT3D_FAILURE;
        }

        // Copy each frame into the mask
        for (auto it = ref_instances.begin(); it != ref_instances.end(); 
                ++it) {

                // Get the SOPInstanceUID
                OFString ref_uid;
                status = (*it)->getReferencedSOPInstanceUID(ref_uid);
                if (status.bad()) {
                        SIFT3D_ERR("read_dso: Failed to get "
                                "ReferencedSOPInstanceUID (%s) \n",
                                status.text());
                        im_free(&maskData);
                        return SIFT3D_FAILURE;
                }
                const char *ref_uid_str = ref_uid.c_str();

                // Check for a matching UID in the images
                auto match = std::find_if(
                        images.begin(), 
                        images.end(), 
                        [ref_uid_str] (Dicom &image) {
                                return image.eqInstance(ref_uid_str);
                        }
                );
                if (match == images.end()) {
                        SIFT3D_ERR("read_dso: No image found with "
                                "SOPInstanceUID %s \n", ref_uid_str);
                        im_free(&maskData);
                        return SIFT3D_FAILURE;
                }

                // Get the z offsets
                const int im_frame_num = std::distance(images.begin(), match);
                const int dso_frame_num = std::distance(ref_instances.begin(), 
                                                        it);

                // Copy the frame into the full mask volume
                if (write_subvolume(*match, mask, im_frame_num, &maskData, 
                        dso_frame_num, dso_frame_num)) {
                        im_free(&maskData);
                        return SIFT3D_FAILURE;
                }
        }

        // Clean up
        im_free(&maskData);
        return SIFT3D_SUCCESS;
}

/* Helper function to read the metadata from a directory of DICOM files */
static int read_dcm_dir_meta(const char *path, std::vector<Dicom> &dicoms) {

        struct stat st;
        DIR *dir;
        struct dirent *ent;

        // Verify that the directory exists
	if (stat(path, &st)) {
                SIFT3D_ERR("read_dcm_dir_cpp: cannot find file %s \n", path);
                return SIFT3D_FAILURE;
	} else if (!S_ISDIR(st.st_mode)) {
                SIFT3D_ERR("read_dcm_dir_cpp: file %s is not a directory \n",
                        path);
                return SIFT3D_FAILURE;
	}

        // Open the directory
        if ((dir = opendir(path)) == NULL) {
                SIFT3D_ERR("read_dcm_dir_cpp: unexpected error opening "
                        "directory %s \n", path);
                return SIFT3D_FAILURE;
        }

        // Get all of the .dcm files in the directory, ignoring DSOs
        dicoms.clear();
        while ((ent = readdir(dir)) != NULL) {

                // Form the full file path
                std::string fullfile(std::string(path) + sepStr + ent->d_name);

                // Check if it is a DICOM file 
                if (im_get_format(fullfile.c_str()) != DICOM)
                        continue;

                // Read the file
                Dicom dicom(fullfile.c_str(), 0);
                if (!dicom.isValid()) {
                        closedir(dir);
                        return SIFT3D_FAILURE;
                }

                // Ignore DSOs
                if (dicom.getType() == DSO)
                        continue;

                // Add the file to the list
                dicoms.push_back(dicom);
        }

        // Release the directory
        closedir(dir);
        
        // Verify that dicom files were found
        if (dicoms.size() == 0) {
                SIFT3D_ERR("read_dcm_dir_cpp: no DICOM files found in %s \n",
                        path);
                return SIFT3D_FAILURE;
        }

        // Sort the slices by their coordinates
        std::sort(dicoms.begin(), dicoms.end()); 

        return SIFT3D_SUCCESS;
}

/* Resize an image to fit a DICOM series. */
static int dcm_resize_im(const std::vector<Dicom> &dicoms, Image *const im) {

        int i, j;

        // Verify that the files are from the same series, and other metadata
        const int num_files = dicoms.size();
        const Dicom &first = dicoms[0];
        const int sortAxis = first.getSortAxis();
        for (int i = 1; i < num_files; i++) {

                const Dicom &dicom = dicoms[i];

                // Verify the SOPSeriesUID
                if (!first.eqSeries(dicom)) {
                        SIFT3D_ERR("read_dcm_dir_cpp: file %s is from a "
                                "different series than file %s \n", 
                                dicom.name(), first.name());
                        return SIFT3D_FAILURE;
                }

                // Verify the sorting axis
                if (dicom.getSortAxis() != sortAxis) {
                        SIFT3D_ERR("read_dcm_dir_cpp: file %s (%d) is sorted "
                                "by a different axis than file %s (%d) \n",
                                dicom.name(), dicom.getSortAxis(), first.name(),
                                sortAxis);
                        return SIFT3D_INCONSISTENT_AXES;
                }
        }

        // Initialize the units to those of the first image
        im->ux = first.getUx(); 
        im->uy = first.getUy();
        im->uz = first.getUz();

        // Compute the spacing between slices
        if (num_files > 1) {

                // Tolerance for spacing discrepancies
                const double tol_spacing = 5E-2;

                // If there are multiple slices, compute the slice spacing by 
                // subtracting the first two images
                const double first_spacing = fabs(first.getSortCoord() - 
                        dicoms[1].getSortCoord());

                // Ensure all the other spacings agree, else throw an error
                for (int i = 0; i < num_files - 1; i++) {

                        const Dicom &prev = dicoms[i];
                        const Dicom &next = dicoms[i + 1];

                        const double spacing = fabs(prev.getSortCoord() - 
                                next.getSortCoord());

                        // Check for duplicates
                        if (spacing == 0.0) {
                                SIFT3D_ERR("read_dcm_dir_cpp: files %s and %s "
                                "have duplicate coordinates in the sorting "
                                "dimension (%d) \n", prev.name(), next.name(),
                                prev.getSortAxis());
                                return SIFT3D_DUPLICATE_SLICES;
                        }

                        // Check the spacing
                        if (fabs(spacing - first_spacing) <= tol_spacing)
                                continue;

                        SIFT3D_ERR("read_dcm_dir_cpp: files %s and %s are "
                                "spaced %fmm apart, which does not follow the "
                                "series spacing of %fmm \n", 
                                prev.name(), next.name(), spacing, 
                                first_spacing);
                        return SIFT3D_UNEVEN_SPACING;
                }

                // Ensure all the slice thicknesses agree, else print a warning
                for (int i = 1; i < num_files; i++) {

                        const Dicom &dicom = dicoms[i];
                        const double thickness = dicom.getSortUnit();

                        if (fabs(first_spacing - thickness) <= tol_spacing)
                                continue; 

                        SIFT3D_ERR("read_dcm_dir_cpp: WARNING--file %s has "
                                "a slice thickness of %fmm, which does not "
                                "agree with the series slice spacing of %fmm\n",
                                dicom.name(), thickness, first_spacing);

                        // No need for multiple warnings
                        break;
                }
               
                SIFT3D_IM_GET_UNITS(im)[sortAxis] = first_spacing;
        }

        // Initialize the output dimensions to those of the first image
        int dims[] = {first.getNx(), first.getNy(), first.getNz()};
        int nc = first.getNc();

        // Verify the dimensions of the other files, counting the total
        // number of slices
        int nSlice = 0;
        for (i = 0; i < num_files; i++) {

                // Get a slice
                const Dicom &dicom = dicoms[i];        

                // Verify the channels
                if (dicom.getNc() != nc) {
                        SIFT3D_ERR("read_dcm_dir_cpp: slice %s "
                                "(%d) does not match the "
                                "channels of slice %s (%d) \n",
                                dicom.name(), dicom.getNc(), 
                                first.name(), nc);
                }

                // Verify the non-sorting dimensions
                for (j = 0; j < 2; j++) {
                        const int firstDim = dims[j];
                        const int sliceDim = dicom.getDim(j);

                        if (sliceDim == firstDim)
                                continue;

                        SIFT3D_ERR("read_dcm_dir_cpp: slice %s "
                                "dimension %d (%d) does not match the "
                                "dimensions of slice %s (%d) \n",
                                dicom.name(), sliceDim,  sortAxis, first.name(),
                                firstDim);
                        return SIFT3D_FAILURE;
                }

                // Count the number of slices
                nSlice += dicom.getDim(sortAxis);
        }

        // Set the dimension of the sorting axis
        dims[sortAxis] = nSlice;

        // Resize the output
        memcpy(SIFT3D_IM_GET_DIMS(im), dims, sizeof(dims));
        im->nc = nc;
        im_default_stride(im);
        if (im_resize(im))
                return SIFT3D_FAILURE;

        return SIFT3D_SUCCESS;
}

/* Copy a Dicom sub-volume into a larger volume. start and end are indices into
 * sub, along the sorting dimension.  */
static int write_subvolume(const Dicom &dicom, Image *const main, 
        const int mainOffset, Image *const sub, const int start, 
        const int end) {

        int mainOffsets[] = {0, 0, 0};
        int starts[] = {0, 0, 0};
        int ends[] = {sub->nx - 1, sub->ny - 1, sub->nz - 1};
        int x, y, z, c, k;
        
        const int sortAxis = dicom.getSortAxis();
        const int mainSortDim = SIFT3D_IM_GET_DIMS(main)[sortAxis];
        const int subSortDim = SIFT3D_IM_GET_DIMS(sub)[sortAxis];

        // Verify the image axis dimensions
        for (k = 0; k < 2; k++) {
                const int axisIdx = dicom.getAxes(k);
                const int mainDim = SIFT3D_IM_GET_DIMS(main)[axisIdx];
                const int subDim = SIFT3D_IM_GET_DIMS(sub)[axisIdx];
                if (mainDim != subDim) {
                        SIFT3D_ERR("write_subvolume: axis %d does not match: "
                                "%d (main) vs %d (sub)", axisIdx, mainDim, 
                                subDim);
                        return SIFT3D_FAILURE;
                }
        }

        // Verify the sorting axis dimensions
        if (start < 0) {
                SIFT3D_ERR("write_subvolume: invalid start %d\n", start);
                return SIFT3D_FAILURE;
        }
        if (end > subSortDim) {
                SIFT3D_ERR("write_subvolume: invalid end %d (maximum %d)\n", 
                        end, subSortDim);
                return SIFT3D_FAILURE;
        }
        if (mainSortDim < mainOffset + end - start) {
                SIFT3D_ERR("write_subvolume: subvolume range (%d, %d) does "
                        "not fit into main volume dimension %d, offset %d, "
                        "axis %d\n", start,  end, mainSortDim, mainOffset, 
                        sortAxis);
                return SIFT3D_FAILURE;
        }

        // Set the starting and ending indices for the sort axis
        starts[sortAxis] = start;
        ends[sortAxis] = end;
        mainOffsets[sortAxis] = mainOffset;

        // Copy the data
        SIFT3D_IM_LOOP_LIMITED_START_C(sub, x, y, z, c, starts[0], ends[0],
                starts[1], ends[1], starts[2], ends[2])
                SIFT3D_IM_GET_VOX(main, x + mainOffsets[0], y + mainOffsets[1],
                        z + mainOffsets[2], c) = 
                        SIFT3D_IM_GET_VOX(sub, x, y, z, c);
        SIFT3D_IM_LOOP_END_C

        return SIFT3D_SUCCESS;
}

/* Helper function to read a directory of DICOM files using C++ */
static int read_dcm_dir_cpp(const char *path, Image *const im) {

        int ret, i, j, sliceCount;

        // Read the DICOM metadata
        std::vector<Dicom> dicoms;
        if (read_dcm_dir_meta(path, dicoms))
                return SIFT3D_FAILURE;

        // Initialize the image volume
        if (ret = dcm_resize_im(dicoms, im))
                return ret;

        // Allocate a temporary image for the slices
        Image slice;
        init_im(&slice);

        // Read the image data
        sliceCount = 0;
        const int num_files = dicoms.size();
        for (i = 0; i < num_files; i++) {

                Dicom dicom = dicoms[i];

                // Read the slice data and write it into the volume
                const int numSlices = dicom.getDim(dicom.getSortAxis());
                if (read_dcm_img(dicom, &slice) ||
                        write_subvolume(dicom, im, sliceCount, &slice, 0, 
                                numSlices - 1)) { 
                        im_free(&slice);
                        return SIFT3D_FAILURE;
                }

                sliceCount += numSlices;
        }
        im_free(&slice);

        return SIFT3D_SUCCESS;
} 

/* Helper function to set meta_new to default values if meta is NULL,
 * otherwise copy meta to meta_new */
static void set_meta_defaults(const Dcm_meta *const meta, 
        Dcm_meta *const meta_new) {
        if (meta == NULL) {
                default_Dcm_meta(meta_new);
        } else {
                *meta_new = *meta;        
        }
}

/* Helper function to write a DICOM file using C++ */
static int write_dcm_cpp(const char *path, const Image *const im,
        const Dcm_meta *const meta, const float max_val) {

#define BUF_LEN 1024
        char buf[BUF_LEN];

        // Ensure the image is monochromatic
        if (im->nc != 1) {
                SIFT3D_ERR("write_dcm_cpp: image %s has %d channels. "
                        "Currently only single-channel images are supported.\n",
                         path, im->nc);
                return SIFT3D_FAILURE;
        }

        // If no metadata was provided, initialize default metadata
        Dcm_meta meta_new;
        set_meta_defaults(meta, &meta_new);

        // Create a new fileformat object
        DcmFileFormat fileFormat;

        // Set the file type to derived
        DcmDataset *const dataset = fileFormat.getDataset();
        OFCondition status = dataset->putAndInsertString(DCM_ImageType, 
                                                         "DERIVED");
        if (status.bad()) {
                std::cerr << "write_dcm_cpp: Failed to set the image type" <<
                        std::endl;
                return SIFT3D_FAILURE;
        }

        // Set the class UID
        dataset->putAndInsertString(DCM_SOPClassUID, 
                UID_CTImageStorage);
        if (status.bad()) {
                SIFT3D_ERR("write_dcm_cpp: Failed to set the SOPClassUID\n");
                return SIFT3D_FAILURE;
        }

        // Set the photometric interpretation
        const char *photoInterp;
        if (im->nc == 1) {
                photoInterp = "MONOCHROME2";
        } else if (im->nc == 3) {
                photoInterp = "RGB";
        } else {
                SIFT3D_ERR("write_dcm_cpp: failed to determine the "
                        "photometric representation for %d channels \n", 
                        im->nc);
                return SIFT3D_FAILURE;
        }
        dataset->putAndInsertString(DCM_PhotometricInterpretation,
                photoInterp);
        if (status.bad()) {
                SIFT3D_ERR("write_dcm_cpp: Failed to set the photometric "
                        "interpretation \n");
                return SIFT3D_FAILURE;
        }

        // Set the pixel representation to unsigned
        dataset->putAndInsertUint16(DCM_PixelRepresentation, 0);
        if (status.bad()) {
                SIFT3D_ERR("write_dcm_cpp: Failed to set the pixel "
                        "representation \n");
                return SIFT3D_FAILURE;
        }

        // Set the number of channels (Samples Per Pixel) and set the planar
        // configuration to interlaced pixels
        assert(SIFT3D_IM_GET_IDX(im, 0, 0, 0, 1) == 
                SIFT3D_IM_GET_IDX(im, 0, 0, 0, 0) + 1);
        snprintf(buf, BUF_LEN, "%d", im->nc);
        dataset->putAndInsertString(DCM_SamplesPerPixel, buf);
        dataset->putAndInsertString(DCM_PlanarConfiguration, "0");

        // Set the bits allocated and stored, in big endian format 
        const unsigned int dcm_high_bit = dcm_bit_width - 1;
        dataset->putAndInsertUint16(DCM_BitsAllocated, dcm_bit_width);
        dataset->putAndInsertUint16(DCM_BitsStored, dcm_bit_width);
        dataset->putAndInsertUint16(DCM_HighBit, dcm_high_bit);
        if (status.bad()) {
                SIFT3D_ERR("write_dcm_cpp: Failed to set the bit widths \n");
                return SIFT3D_FAILURE;
        }

        // Set the patient name
        status = dataset->putAndInsertString(DCM_PatientName, 
                meta_new.patient_name);
        if (status.bad()) {
                SIFT3D_ERR("write_dcm_cpp: Failed to set the patient name\n");
                return SIFT3D_FAILURE;
        }

        // Set the patient ID
        status = dataset->putAndInsertString(DCM_PatientID,
                meta_new.patient_id);
        if (status.bad()) {
                SIFT3D_ERR("write_dcm_cpp: Failed to set the patient ID \n");
                return SIFT3D_FAILURE;
        }

        // Set the study UID
        status = dataset->putAndInsertString(DCM_StudyInstanceUID,
                meta_new.study_uid);
        if (status.bad()) {
                SIFT3D_ERR("write_dcm_cpp: Failed to set the "
                        "StudyInstanceUID \n");
                return SIFT3D_FAILURE;
        }

        // Set the series UID
        status = dataset->putAndInsertString(DCM_SeriesInstanceUID,
                meta_new.series_uid);
        if (status.bad()) {
                SIFT3D_ERR("write_dcm_cpp: Failed to set the "
                        "SeriesInstanceUID \n");
                return SIFT3D_FAILURE;
        }

        // Set the series description
        status = dataset->putAndInsertString(DCM_SeriesDescription,
                meta_new.series_descrip);
        if (status.bad()) {
                SIFT3D_ERR("write_dcm_cpp: Failed to set the series "
                        "description \n");
                return SIFT3D_FAILURE;
        }

        // Set the instance UID
        status = dataset->putAndInsertString(DCM_SOPInstanceUID, 
                meta_new.instance_uid);
        if (status.bad()) {
                SIFT3D_ERR("write_dcm_cpp: failed to set the "
                        "SOPInstanceUID \n");
                return SIFT3D_FAILURE;
        }

        // Set the dimensions
        OFCondition xstatus = dataset->putAndInsertUint16(DCM_Rows, im->ny); 
        OFCondition ystatus = dataset->putAndInsertUint16(DCM_Columns, im->nx);
        snprintf(buf, BUF_LEN, "%d", im->nz);
        OFCondition zstatus = dataset->putAndInsertString(DCM_NumberOfFrames,
                buf);
        if (xstatus.bad() || ystatus.bad() || zstatus.bad()) {
                SIFT3D_ERR("write_dcm_cpp: Failed to set the dimensions \n");
                return SIFT3D_FAILURE;
        }

        // Set the instance number
        snprintf(buf, BUF_LEN, "%u", meta_new.instance_num);
        status = dataset->putAndInsertString(DCM_InstanceNumber, buf);
        if (status.bad()) {
                SIFT3D_ERR("write_dcm_cpp: Failed to set the instance "
                        "number \n");
                return SIFT3D_FAILURE;
        }

        // Set the ImagePositionPatient vector
        const double imPosX = static_cast<double>(im->nx - 1) * im->ux;
        const double imPosY = static_cast<double>(im->ny - 1) * im->uy;
        const double imPosZ = static_cast<double>(meta_new.instance_num) * 
                              im->uz;
        snprintf(buf, BUF_LEN, imPosPatientFmt, imPosX, imPosY, imPosZ);
        status = dataset->putAndInsertString(DCM_ImagePositionPatient, buf);
        if (status.bad()) {
                SIFT3D_ERR("write_dcm_cpp: Failed to set the "
                        "ImagePositionPatient vector \n");
                return SIFT3D_FAILURE;
        }

        // Set the ImageOrientationPatient vectors
        const float imageOriPatient[] = {1., 0, 0, 0, 1., 0};
        snprintf(buf, BUF_LEN, imOriPatientFmt, 
                imageOriPatient[0], 
                imageOriPatient[1],
                imageOriPatient[2],
                imageOriPatient[3],
                imageOriPatient[4],
                imageOriPatient[5]);
        status = dataset->putAndInsertString(DCM_ImageOrientationPatient, buf);
        if (status.bad()) {
                SIFT3D_ERR("write_dcm_cpp: Failed to set the "
                        "ImageOrientationPatient vectors \n");
                return SIFT3D_FAILURE;
        }

        // Set the slice location
        snprintf(buf, BUF_LEN, "%f", imPosZ);
        status = dataset->putAndInsertString(DCM_SliceLocation, buf);
        if (status.bad()) {
                SIFT3D_ERR("write_dcm_cpp: Failed to set the slice "
                        "location \n");
                return SIFT3D_FAILURE;
        }

        // Set the pixel spacing
        snprintf(buf, BUF_LEN, pixelSpacingFmt, im->ux, im->uy);
        status = dataset->putAndInsertString(DCM_PixelSpacing, buf);
        if (status.bad()) {
                SIFT3D_ERR("write_dcm_cpp: Failed to set the pixel "
                        "spacing \n");
                return SIFT3D_FAILURE;
        }

        // Set the aspect ratio
        snprintf(buf, BUF_LEN, "%f\\%f", im->ux, im->uy);
        status = dataset->putAndInsertString(DCM_PixelAspectRatio, buf);
        if (status.bad()) {
                SIFT3D_ERR("write_dcm_cpp: Failed to set the pixel aspect "
                        "aspect ratio \n");
                return SIFT3D_FAILURE;
        }

        // Set the slice thickness
        snprintf(buf, BUF_LEN, "%f", im->uz);
        status = dataset->putAndInsertString(DCM_SliceThickness, buf);
        if (status.bad()) {
                SIFT3D_ERR("write_dcm_cpp: Failed to set the slice "
                                "thickness \n");
                return SIFT3D_FAILURE;
        }

        // Count the number of pixels in the image
        unsigned long numPixels = SIFT3D_IM_GET_DIMS(im)[0];
        for (int i = 1; i < IM_NDIMS; i++) {
                numPixels *= SIFT3D_IM_GET_DIMS(im)[i];
        }

        // Get the image scaling factor
        const float dcm_max_val = static_cast<float>(1 << dcm_bit_width) - 1.0f;
        const float im_max = max_val < 0.0f ? im_max_abs(im) : max_val;
        const float scale = im_max == 0.0f ? 1.0f : dcm_max_val / im_max;

        // Render the data to an 8-bit unsigned integer array
        assert(dcm_bit_width == 8);
        assert(fabsf(dcm_max_val - 255.0f) < FLT_EPSILON);
        uint8_t *pixelData = new uint8_t[numPixels];
        int x, y, z, c;
        SIFT3D_IM_LOOP_START_C(im, x, y, z, c)

                const float vox = SIFT3D_IM_GET_VOX(im, x, y, z, c);

                if (vox < 0.0f) {
                        SIFT3D_ERR("write_dcm_cpp: Image cannot be "
                                "negative \n");
                        return SIFT3D_FAILURE;
                }

                pixelData[c + x + y * im->nx + z * im->nx * im->ny] =
                        static_cast<uint8_t>(vox * scale);
        SIFT3D_IM_LOOP_END_C

        // Write the data
        status = dataset->putAndInsertUint8Array(DCM_PixelData, pixelData, 
                numPixels);
        delete[] pixelData;
        if (status.bad()) {
                SIFT3D_ERR("write_dcm_cpp: failed to set the pixel data \n");
                return SIFT3D_FAILURE;
        }

        // Choose the encoding format
#if 0
        DJEncoderRegistration::registerCodecs();
        const E_TransferSyntax xfer = EXS_JPEGProcess14SV1TransferSyntax;
        DJ_RPLossless rp_lossless;
        status = dataset->chooseRepresentation(xfer, &rp_lossless);
#else
        const E_TransferSyntax xfer = EXS_LittleEndianExplicit;
        dataset->chooseRepresentation(xfer, NULL);
#endif
        if (!dataset->canWriteXfer(xfer)) {
                SIFT3D_ERR("write_dcm_cpp: Failed to choose the encoding "
                        "format \n");
                return SIFT3D_FAILURE;
        }

        // Force the media storage UIDs to be re-generated by removing them
        dataset->remove(DCM_MediaStorageSOPClassUID);
        dataset->remove(DCM_MediaStorageSOPInstanceUID);

        // Save the file
        status = fileFormat.saveFile(path, xfer);
        if (status.bad()) {
                SIFT3D_ERR("write_dcm_cpp: failed to write file %s (%s) \n",
                        path, status.text());
                return SIFT3D_FAILURE;
        }

        return SIFT3D_SUCCESS;
#undef BUF_LEN
}

/* Helper function to write an image to a directory of DICOM files using C++ */
static int write_dcm_dir_cpp(const char *path, const Image *const im,
        const Dcm_meta *const meta) {

        Image slice;

        // Initialize C intermediates
        init_im(&slice);

        // Initialize the metadata to defaults, if it is null 
        Dcm_meta meta_new;
        set_meta_defaults(meta, &meta_new);

        // Get the number of leading zeros for the file names
        const int num_slices = im->nz;
        const int num_zeros = static_cast<int>(ceil(log10(
                static_cast<double>(num_slices))));

        // Form the printf format string for file names
#define BUF_LEN 16
        char format[BUF_LEN];
        snprintf(format, BUF_LEN, "%%0%dd.%s", num_zeros, ext_dcm); 
#undef BUF_LEN

        // Resize the slice buffer
        slice.nx = im->nx; 
        slice.ny = im->ny;
        slice.nz = 1;
        slice.nc = im->nc;
        im_default_stride(&slice);
        if (im_resize(&slice)) {
                im_free(&slice);
                return SIFT3D_FAILURE;
        }

        // Copy the units to the slice
        memcpy(SIFT3D_IM_GET_UNITS(&slice), SIFT3D_IM_GET_UNITS(im), 
                IM_NDIMS * sizeof(double));

        // Get the maximum absolute value of the whole image volume
        const float max_val = im_max_abs(im);

        // Write each slice
        for (int i = 0; i < num_slices; i++) {

                // Form the slice file name
#define BUF_LEN 1024
                char buf[BUF_LEN];
                snprintf(buf, BUF_LEN, format, i);

                // Form the full file path
                std::string fullfile(path + sepStr + buf);

                // Copy the data to the slice
                int x, y, z, c;
                SIFT3D_IM_LOOP_START_C(&slice, x, y, z, c)
                        SIFT3D_IM_GET_VOX(&slice, x, y, z, c) =
                                SIFT3D_IM_GET_VOX(im, x, y, i, c);
                SIFT3D_IM_LOOP_END_C

                // Generate a new SOPInstanceUID
                dcmGenerateUniqueIdentifier(meta_new.instance_uid, 
                        SITE_INSTANCE_UID_ROOT); 

                // Set the instance number
                const unsigned int instance = static_cast<unsigned int>(i + 1);
                meta_new.instance_num = instance;

                // Write the slice to a file
                if (write_dcm(fullfile.c_str(), &slice, &meta_new, max_val)) {
                        im_free(&slice);
                        return SIFT3D_FAILURE;
                }
        }

        // Clean up
        im_free (&slice);

        return SIFT3D_SUCCESS;
}

#endif
